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as predicted by Lagrangian perturbation approximations up to the third order. 
The first-order approximation is equivalent to the "Zel'dovich approximation" 
for the type of initial data analyzed. The study is performed for two simple 
models which allow studying of typical features of the clustering process in 
the early non-linear regime. We calculate the initial perturbation potentials 
as solutions of Poisson equations algebraically, and automate this calculation 
for a given initial random density field. The presented models may also be 
useful for other questions addressed to Lagrangian perturbation solutions and 
for the comparison of different approximation schemes. In an accompanying 
paper we investigate a detailed comparison with various N-body integrators 
using these models (Karakatsanis & Buchert 1995). 

Results of the present paper include the following: 1. The collapse is acceler- 
ated significantly by the higher-order corrections confirming previous results 
by Moutarde et al. (1991); 2. the spatial structure of the density patterns pre- 
dicted by the "Zel'dovich approximation" differs much from those predicted by 
the second- and third-order Lagrangian approximations; 3. Second-order ef- 
fects amount to internal substructures such as "second generation" -pancakes, 
-filaments and -clusters, as are also observed in N-body simulations; 4. The 
third-order effect gives rise to substructuring of the secondary mass-shells. 
The hierarchy of shell-crossing singularities that form features small high- 
density clumps at the intersections of caustics which we interprete as gravita- 
tional fragmentation. 



at high resolution 





I We present high-spatial resolution studies of the density field 
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1. Introduction 



It is commonly appreciated that Lagrangian perturbation solutions provide useful models of 
large-scale structure. Comparison with numerical simulations have put them into a strong 
position in the list of currently discussed analytical or semi-analytical approximations (see: 
Melott 1994 for a summary). Lagrangian perturbation schemes have been optimized by 
smoothing the high-frequency end of the power spectrum of density inhomogeneities such 
that they are capable of replacing N-body integrators above some scale close, but smaller 
than the non-linearity scale (i.e., where the r.m.s. density contrast is of order unity) (Coles 
et al. 1993, Melott eA al. 1994, 1995, Bouchet eA al. 1995, Sathyaprakhash eA al. 1995, 
Weifi eA al. 1995). While their application to pancake models, i.e., models with a large 
high-frequency cutoff, demonstrates an excellent performance up to the epoch when shell- 
crossing singularities in the cosmic flow develop (Buchert eA al. 1994), their application 
to later non-linear stages fails unless the initial data are smoothed to avoid substantial 
post-singularity evolution. This way the large-scale structure is restored, and small-scale 
features arise due to the collapse of the waves which were left in the initial data. 

In the present work we want to look in more detail at the collapsing structures around 
the epoch of shell-crossing on smaller scales by using high-resolution techniques described 
by Buchert & Bartelmann (1991) (however, here, we do not interpolate initial data). The 
present study can be viewed in line with previous high-resolution studies of pancakes 
(Buchert 1989a,b, Melott & Shandarin 1990, Beacom et al. 1991 (2D), and Buchert & 
Bartelmann 1991 (2D and 3D)). 

The applicability of the Lagrangian approximations has been tested in previous work 
on the basis of cross-correlation statistics of density flelds in which the internal substruc- 
tures are not resolved. We here address the question which substructures are predicted by 
these approximations and we shall single out flrst-, second-, and third-order effects in the 
evolution of caustics in the density fleld. The work by Alimi et al. (1990), Moutarde et 
al. (1991) and the comparison of Lagrangian perturbation solutions with the spherically 
symmetric solution by Munshi et al. (1994) comprise steps in this direction. 

We have taken care of the precision with which we realize the Lagrangian schemes. 
Thus far, these analytical models have to be realized numerically to set up the initial 
data in Fourier space. In particular, the third-order model provides a complication, since 
products of derivatives of perturbation potentials at flrst and second order (as solutions 
of Poisson equations) form the input for the third-order perturbation potentials, which 
describe interaction of perturbations (see the next subsection for details). It is therefore 
desirable to control this realization of initial data in an optimal way to minimize numerical 
uncertainties. We did this by calculating the perturbations fully analytically. We have 
also automated the process of flnding the perturbation potentials from a single given 



2 



initial velocity potential, or density field, respectively, by using algebraic manipulation 
systems. This procedure is suitable for spectra with not too much modes like in models 
with coherence length. Another advantage of this analytical procedure is the possibility 
of improving particle number, since we are not limited by storage as in the case of FFT 
realizations. We therefore can present realizations using 1024^ particles. The structures 
shown can only be seen at resolutions higher than 256^ particles. 

We start with the derivation of a simple plane-wave model attempted earlier by 
Moutarde et al. (1991) (see also Alimi eA al. 1990). Their third-order solution was 
derived just for this model. (However, this solution did not pass a test we did by inserting 
it into the Euler-Poisson system in Lagrangian form.) We, here, proceed differently. We 
start from the generic third-order solution given by Buchert (1994) and insert the plane- 
wave model as a special case. Since both the generic model and the special model have 
been checked to solve the original equations (by using algebraic manipulation systems), we 
are confident in our calculations. Besides automating algebraically the derivation of the 
potentials as mentioned above, we also exemplify the use of a set of local forms given by 
Buchert & Ehlers (1993) and Buchert (1994) in the case of the Moutarde eA al. problem 
(see the APPENDIX). We stick to that model first, since it is simple and already shows 
the principal features of the gravitational collapse we are interested in. Also in other work 
on related subjects this model is useful as an example (Mo & Buchert 1990, Matarrese eA 
al. 1992), and can be used as a toy-model to compare different approximation schemes. 
We then move to generic initial data, i.e., data with no symmetry, but restricted to a small 
enough box to assure the resolution of patterns we are interested in. 
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2. A generic third-order solution 



Let us recall the class of third-order solutions on which we base our models. We require 
that, initially, the peculiar-velocity u(X^t) to be proportional to the peculiar-acceleration 

u(X,to) = w(X,to)to , (1) 

where we have defined the fields as usual (compare Peebles 1980, Buchert 1992). Hence- 
forth, we denote the peculiar-velocity potential at the initial time to by iS, to) =: Vo<S, 
and the peculiar-gravitational potential at to by ^, w(X^to) =: — Vo^, where Vo denotes 
the nabla operator with respect to the Lagrangian coordinates X. The restriction (1) 
has proved to be appropriate for the purpose of modeling large-scale structure, since the 
peculiar-velocity field tends to be parallel to the gravitational peculiar-field strength after 
some time, related to the existence of growing and decaying solutions in the linear regime 
(Buchert & Ehlers 1993, Buchert 1994). 

With a superposition ansatz for Lagrangian perturbations of an Einstein-de Sitter 
background the following mapping q = F(X^a) as irrotational solution of the Euler- 
Poisson system in Lagrangian form up to the third order in the perturbations from ho- 
mogeneity has been obtained (Buchert 1994). F defines the displacement map from La- 
grangian coordinates X to Eulerian coordinates q which are comoving with the unperturbed 
Hubble-fiow; the general set of initial conditions (^(X),iS(X)) is restricted according to 
S = —(pto (see equation (1)); a(t) = (t/to)^'^^: 

F = X + q,{a)VoS'^'\X) + q2{a)VoS'^^\X) 
+ qlia) VoS^'^'HX) + qlia) VoS^'^HX) - ql{a) Vo X ^(^^^(X) , (2) 

with: 
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where the initial displacement vectors have to be constructed by solving seven elliptic 
boundary value problems (summation over repeated indices; i,j,k = 1,2,3 with cyclic or- 
dering). 

AoS^'^ = I{S,,,k)to , (2/) 
Ao<S(^) = 2//(4^i) , (2^,) 



a,b,c 



An important remark relevant to any realization of the solution (2) concerns the possi- 
bility of setting S^^^ = Sto without loss of generality, if initial data are spatially periodic 
(compare Buchert 1995b, Ehlers & Buchert 1995 for details and proofs). With this setting, 
the first-order solution reduces to the well-known "Zel'dovich approximation" (Zel'dovich 
1970, 1973; Buchert 1992), which then assumes its familiar local form. Also, the trun- 
cated third-order model (i.e., neglecting interaction terms) is then, although of course 
non-locally, expressible in terms of the initial data (compare eqs. (2f-h)). 

The scalar potential iS*^^^-* and the vector potential iS*^^^-* generate interaction among 
the first- and second-order perturbations. The general interaction term is not purely lon- 
gitudinal: inspite of the irrotationality of the flow in Eulerian space, vorticity is generated 
in Lagrangian space starting at the third order for this set of intial data. For more general 
initial data, this happens already at second order. As our analysis of the solution will 
show, it has sense to include the interaction term iS*^^^-* only, neglecting the transverse part 
altogether. However, as will be demonstrated, keeping only the generating function iS*^^"-* 
results in a density pattern, which is not an adequate generalization of the second-order 
approximation. This "truncated third-order" model has been proposed in (Buchert 1994) 
as the "main body" of the perturbation sequence in the early nonlinear regime, since all 
higher-order solutions are made up of interaction terms among the perturbation poten- 
tials. A closer look at the features presented in this work shows that the third-order model 
has to he run with the interaction term S^'^^\ 
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3. Special clustering models in closed form 



In the APPENDIX we demonstrate how to construct special models by using "local forms" 
for the displacement vectors. Although analytically interesting, this procedure is cumber- 
some if applied to more complex initial data. In this section we describe how we can 
automate the process of finding closed form expressions for the perturbation potentials. 

In general, we are interested in a class of initial data which can be represented by a finite 
Fourier sum of plane waves having random amplitudes and random phases. The random 
variables can, e.g., be specified in terms of a power spectrum of a Gaussian random density 
field. Usually, such initial conditions are generated by FFT (Fast Fourier Transform), a 
method which was also used to realize the generic model in (Buchert et al. 1994, Melott eA 
al. 1995, Weifi eA al. 1995). However, there are two limitations of this method which both 
restrict the power of spatial resolution, an advantage which is in principle offered by ana- 
lytical solutions. One of these limitations is due to the limited CPU storage for employing 
the FFT routine, the other is due to a lack of precision which may arise by constructing 
the initially small displacements from a given density field, or by interpolating the particle 
displacements into a smooth density field (using, e.g., CIC binning), respectively. As an 
alternative, we suggest to solve the Poisson equations in eqs. (2) algebraically by compar- 
ing the coefficients of Fourier sums in the source terms and the perturbation potentials. 
This way the solutions can be calculated to high accuracy without hitting on CPU storage 
limitations. Since the model is a one-timestep mapping, the CPU time needed for the 
realization is still reasonably small (512^ particles require CPU times of a few hours for 
the generic model discussed below). However, we admit that the algebraic procedure to 
solve for all seven perturbation potentials in (2) is still limited by the CPU storage for the 
algebraic program, and the compilation time of, e.g., plot routines can be large for a large 
number of Fourier modes. Using the manipulation system Mathematical we are easily 
able to construct all perturbation potentials for 50 Fourier modes on a workstation with 
256M storage. The results obtained with this method have also been checked to solve the 
original equations using two algebraic manipulation systems (Reduce and Mathematica). 

For the special models constructed algebraically in this way we have also run the 
previous code (using FFT), which constructs displacements from given density fields (A.G. 
Weifi , priv. comm.), and found as expected that the result is a slightly smoothed variant of 
a direct calculation pursued in the present work. At the same time, this was an independent 
check of the third-order program used in previous work (compare Weifi & Buchert 1993). 

Besides the special model given in the APPENDIX (Model I), i.e., for the initial 
potential 

S\^\X):= --^ (a^ cos{27rX) + ay cos{27rY) + a, cos{27rZ)) , e = | , (3a) 
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we have analyzed a generic model (Model II) with the following initial potential (here, the 
coordinates are normalized by 27r): 

sff{X) := 0.1953 4.82sin(-X-F)+3.95cos(-X-F)+8.82sin(-X-Z)+2.5cos(-X-Z) 

+5.32 sin(-X) + 2.11 cos(-X) + 3.29 sin(-X + Z) + 1.83 cos(-X + Z) + 6.7 sin(-X + Y) 
+4.05 cos(-X + Y) + 6.92 sin(-F - Z) + 1.2 cos(-F - Z) + 3.8 sin(-F) + 4.77 cos(-F) 
+ 1.8 sin(-F + Z) + 4.58 cos(-F + Z) + 1.29 sin(-Z) + 6.6 cos(-Z) + 8.25 sin(Z) 
+4.77 cos(-Z) + 3.67 sin(F - Z) + 3.48 cos(F - Z) + 2.4 sin(F) + 6.02 cos(F) 
+7.86sin(F + Z) + 6.64cos(F + Z) + 9.33sin(X -Y) + 0.87cos(X - Y) 
+0.56 sin(X - Z) + 4.48 cos(X - Z) + 3.4 sin(X) + 5.77 cos(X) + 4.54 sin(X + Z) 

+4.46 cos(X + Z) + 9.8 sin(X + F) + 3.13 cos(X + F) . (36) 

The coefficients have been determined by the requirement that the power spectrum had the 
slope +1 down to the smallest wavelength, and the r.m.s. density contrast had the same 
value as Model I. Model I is the model studied by Moutarde et al. (1991); it has also been 
used by Mo & Buchert (1990) (at first order) as a statistical toy-model, and by Buchert 
& Ehlers (1993) (at second order) to demonstrate secondary shell-crossings; Matarrese 
eA al. (1992) and Kate Croudace (priv. comm.) have compared general relativistic with 
Newtonian dynamics with the help of this model. 

All seven perturbation potentials and the corresponding displacement vectors are listed in 
the APPENDIX for Model I. For Model II the potentials and the displacement vectors can 
be obtained on request. 
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4. High-resolution Studies 



We present high-resolution studies of the density field as predicted by the Lagrangian 
schemes for both models. This is done by collecting 1024^ trajectories into a (comoving) 
Eulerian grid of 512^ cells for Model I (512^ into a grid of 256^ for Model II) (for the 
method see Buchert & Bartelmann 1991). 

Figure 1 displays three evolution stages of the density field predicted by Model I for 
the first-, second-, and third-order perturbation solutions. (Initial data were given at 
Zi = 1000; a.{zi ) = 1.) A manifest feature is the delay of the collapse time for perturbation 
solutions at different orders; higher-order corrections significantly accelerate the collapse. 
This result was already stated by Moutarde eA al. (1991). To compare the spatial patterns 
for the different orders, we can roughly compare the density fields "diagonally" in Fig.l 
(this way of comparison will be discussed quantitatively in a forthcoming paper: Karakat- 
sanis & Buchert 1995): while the first-order solution (the "Zel'dovich approximation") 
carries mainly kinematical information beyond the epoch of shell-crossing, the second- 
order solution modifies the shape of the first mass-shell after crossing and generates a 
second mass-shell as well as secondary sheets and filaments inside the first structures (in 
agreement with the previous study of the trajectory field by Buchert & Ehlers 1993); the 
third-order correction redistributes mass inside the two mass-shells as well as in sheets 
and filaments. 

Figure 2 displays the third-order density field for another color coding which is useful 
to separate the different parts of the third-order corrections in the solution (2): we infer 
that the transverse part of the "interaction term" (2j,k,l) is not of crucial importance 
and might be neglected, it merely deconcentrates the inner mass-shell more from the 
center (which can be seen by comparing full third-order with or without transverse part, 
or "truncated third-order" with or without transverse part). However, to neglect the 
"interaction terms" altogether results in a pattern which is further away from the second- 
order approximation than the full third-order approximation. The outer caustic is even 
absent. This indicates that the third-order approximation without "interaction terms" 
is not useful, the "main body" of the perturbation sequence is not a good model as was 
speculated in (Buchert 1994). 

We continue by looking at Figure 3 which presents the density field of Model II for 
the different orders at a late evolution stage, late, because we then are able to separate 
the different structures visually which appear much earlier in the evolution. Again the 
features quoted above are visible, the collapse is delayed by a huge factor in the first-order 
( "Zel'dovich-" )approximation. Also, the similarity between second- and third-order is 
striking, while the first-order model lacks some internal structures, which can be attributed 
to secondary shell-crossing events (a second-order effect). 
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An interesting aspect of these high-resolution studies relates to a new interpretation of 
the longstanding "fragmentation problem" in classical pancake theory: we appreciate small 
"fragments" sitting at the intersection of caustics (see Figs. 1-3). Since a finite resolution 
brings the density to a finite value, these "fragments" show up as almost spherical blobs 
with potential wells that have about 2 times more height than the potential of the mass- 
shells. In realistic situations, physical processes at the location of caustics and velocity 
dispersion in a dark coUisionless component will do a similar job. We may interprete this 
phenomenon as "gravitational fragmentation": although the initial fluctuation is coherent 
like in pancake models, the collapse process forms fragments on a substantially smaller 
scale. This interpretation is appropriate, if gravity is the dominating interaction related 
to the existence of a mass dominating dark matter component in the Universe. It has far 
reaching consequences in a self-gravitating medium, since we expect the phenomenon of 
multiple shell-crossing (termed "non-dissipative gravitational turbulence by Gurevich & 
Zybin 1988a, b) to continue down to smaller and smaller scales yielding a hierarchy of nested 
caustics. This has been demonstrated in a two-dimensional simulation by Doroshkevich et 
al. (1980). Since further and further mass-shells are generated in the center of a cluster, 
more and more caustics are simultaneously present and consequently create a huge number 
of "fragments" as an internal organization of mass-shells. This way, a cluster naturally 
creates a gravitational potential which is distinctly rippled and thereby prepares the sites 
for galaxy formation: we expect the baryons to preferentially drop into these "fragments". 

Although this consideration has to remain premature at this stage, we think that 
"gravitational fragmentation" as we describe it is a generic effect in gravitational clustering 
and should be taken seriously as soon as a dark matter component dominates the matter 
density. The fact that we need high-spatial resolution studies to uncover these fragments 
explains their absence in the literature. It is interesting to note here that another type of 
fragments appeared in a two-dimensional numerical simulation at high resolution (Melott 
& Shandarin 1990), which results from a redistribution of mass inside fllaments (compare 
their plot with the fllaments in the third-order approximation in Fig.l). 

The detailed study of caustic metamorphoses begun by Arnol'd eA al. (1982) for the 
"Zel'dovich-approximation" in two spatial dimensions will provide the necessary insight 
to further understand this phenomenon. We have continued this study in three spatial 
dimensions (Buchert eA al. 1995a, b); for an overview see (Buchert 1995a). 

Acknowledgments: We would like to thank Jean-Michel Alimi for discussions, and Arno 
Weifi for constructive contributions and discussions. GK and TB were both supported by the "Son- 
derforschungsbereich 375-95 fiir Astro-Teilchenphysik der Deutschen Forschungsgemeinschaft" . 
All authors acknowledge financial support from a DAAD exchange program. 



9 



References 



Alimi J.-M., Bouchet F.R., Pellat R., Sygnet J.-F, Moutarde F. (1990): Ap.J. 354, 3. 

Arnol'd V.I., Shandarin S.F., Zel'dovich Ya.B. (1982): Geophys. Astrophys. Fluid. Dyn. 
20, 111. 

Beacom J.F., Dominik K.G., Melott A.L., Perkins S.P., Shandarin S.F. (1991): Ap.J. 372, 
351. 

Buchert T. (1989a): Astron. Astrophys. 223, 9. 

Buchert T. (1989b): Rev. Mod. Astron. 2, 267. 

Buchert T., Bartelmann M. (1991): Astron. Astrophys. 251, 389. 

Buchert T. (1992): M.N.R.A.S. 254, 729. 

Buchert T., Ehlers J. (1993): M.N.R.A.S. 264, 375. 

Buchert T. (1994): M.N.R.A.S. 267, 811. 

Buchert T. (1995a): in: Formation and Evolution of Galaxies and Structures - The coming 
decade^ eds.: Z.-L. Zou, Y. Chen, P.-W. Ji, Astrophysics Reports, Pub. Beijing Astron. 
Obs. 1, pp. 59-70. 

Buchert T. (1995b): in: Proc. lOP 'Enrico Fermi\ Course CXXXII (Dark Matter in the 
Universe), Varenna 1995, eds.: S. Bonometto, J. Primack, A. Provenzale, in press. 

Buchert T., Melott A.L., Weifi A.G. (1994): Astron. Astrophys. 288, 349. 

Buchert T., Shandarin S.F., Weifi A.G. (1995a): in preparation. 

Buchert T., Shandarin S.F., Weifi A.G. (1995b): in preparation. 

Coles P., Melott A.L., Shandarin S.F. (1993): M.N.R.A.S. 260, 765. 

Doroshkevich A.G., Kotok E.V., Novikov I.D., Polyudov A.N., Shandarin S.F., Sigov Yu.S. 
(1980): M.N.R.A.S. 192, 321. 

Ehlers J., Buchert T. (1995): in preparation. 

Gurevich A.V., Zybin K.P. (1988a): Sov. Phys. JETP 67, 1. 

Gurevich A.V., Zybin K.P. (1988b): Sov. Phys. JETP 67, 1957. 

Karakatsanis G., Buchert T. (1995): Astron. Astrophys., to be submitted. 

Matarrese S., Lucchin F., Moscardini L., Saez D. (1992): M.N.R.A.S. 259, 437. 

Melott A.L. (1994): Ap.J. Lett. 426, Lll. 

Melott A.L., Shandarin S.F. (1990): Nature 346, 633. 

Melott A.L., Pellman T.F., Shandarin S.F. (1994): M.N.R.A.S. 269, 626. 

Melott A.L., Buchert T., Weifi A.G. (1995): Astron. Astrophys. 294, 345. 

Mo H.J., Buchert T. (1990): Astron. Astrophys. 234, 5. 

Moutarde F., Alimi J.-M., Bouchet F.R., Pellat R., Ramani A. (1991): Ap.J. 382, 377. 



10 



Munshi D., Sahni V., Starobinsky A. A. (1994): Ap.J. 436, 517. 

Peebles P.J.E. (1980): The Large Scale Structure of the Universe^ Princeton Univ. Press. 

Sathyaprakash B.S., Sahni V., Munshi D., Pogosyan D., Melott A.L. (1995): M.N.R.A.S., 
in press. 

Weifi A.G., Buchert T. (1993): in: 9th lAP conference: Cosmic Velocity Fields, eds.: F.R. 
Bouchet, M. Lachieze-Rey, Editions Frontieres, pp. 517-519. 

Weifi A.G., Gottl5ber S., Buchert T. (1995): M.N.R.A.S., in press. 

Zel'dovich Ya.B. (1970): Astron. Astrophys. 5, 84. 

Zel'dovich Ya.B. (1973): Astrophysics 6, 164. 

Figure Captions 

Figure 1: Three stages at expansion factors (a = 1000,1200,1500) are shown for the 
first-order approximation (the "Zel'dovich approximation"), (top), the second-order ap- 
proximation (middle), and the third-order approximation (bottom). The initial condition 
is a special periodic function which maps principal elements of the large-scale structure 
such as sheets, filaments and clusters. At a stage shortly after the first shell-crossing (in 
this normalization at a = 1000) the three approximations mainly differ in their predic- 
tion of the collapse time. The higher-order corrections accelerate the collapse significantly 
and generate "second generation" pancakes, filaments and clusters. If we renormalize the 
amplitudes such that the collapse occurs at the same instant in all three approximations, 
then we can read the figures in a diagonal manner, i.e., the second stage in the upper row 
roughly corresponds to the first stage in the middle row, etc. . 

Figure 2: A zoom (1/4 of the box) into the density pattern of the third-order approx- 
imation for Model I is shown with a different color coding. The contributions to the 
third-order effect have been splitted into a part (upper left) which belongs to the "main 
body" (see Section 2), a longitudinal interaction effect added to it (upper right), and a 
transverse interaction effect added to it (lower right). The full third-order approximation 
is shown in the lower left corner. 

Figure 3: A comparison similar to Figure 1, however, for the generic clustering model. 
The density pattern predicted by the first-order ("Zel'dovich-)approximation is shown in 
the upper left panel, that for second-order in the upper right, for third-order in the lower 
left, and for third-order without "interaction terms" in the lower right panel. The same 
effects as quoted for Model I can be seen. 
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APPENDIX 



A.I. The first order displacement vector 

In this section we exemplify how special solutions to (2) can be constructed by using a set of local forms 
which provide first integrals of the seven Poisson equations (2f-l). Let us consider a simple plane-wave model 
for the initial peculiar-velocity potential S; iS*-^-* := Sto, which was studied earlier by Moutarde et al. (1991): 

:= -— ^ (a^cos(27rX) + ttj, cos(27rY) + cos(27rZ)) . (A. la) 

The amplitude e plays the role of the perturbation parameter here and is related to the total amplitude a of 
the density contrast 6 := as a = |£. The amplitudes ax,C(y, ctz allow for triaxial deformations of the 

model; one has to choose aj. -\- a'y -\- aj, = 1 in order to keep the r.m.s. amplitude of b the same. In this paper 
we shall use = l,ay = l,a^ = 1, since different amplitudes give no further information about internal 
structures of the model. Although the model (A. la) is simple, it has no symmetries which destroy the generic 
feature of the singularites formed like plane or spherical symmetry would do. The structure of the cluster 
formed will only retain refiection and rotational symmetries manifest in the potential (A. la) for our choice of 
amplitudes. As, e.g., demonstrated in (Buchert & Ehlers 1993) for a similar two-dimensional model, we have 
with models like (A. la) the possibility of studying principal kinematical features of a generic collapse such 
as the formation of cusped caustics, interconnected network structures, infall of matter onto the cluster. 
Additionally, internal differentiation of a multi-stream system resulting in a hierarchy of shell-crossings, 
which are attributed to a generic feature of a gravitational collapse, can be demonstrated nicely with this 
model. The model has periodic boundary conditions which makes it accessible for numerical treatment. 

From (A. la) we have for the first order displacement vector: 

(a^ sin(27rA') \ 
a,sin(27rY) . {A.lh) 
sin(27r2') / 

We now scetch a procedure how to construct the higher-order potentials from this initial condition. The 
procedure is based on a list of local forms given by Buchert & Ehlers (1993) and Buchert (1994), which 
are, roughly speaking, first integrals of the quadratic and cubic source terms in the Poisson equations of the 
solution (2). These integrals only hold for special classes of initial data, although they might also be useful 
as approximations for generic initial data. For the potential (A. la) it turns out that it belongs to the class 
of initial data which, for all orders, admits such first integrals. 

A. II. Tlie second order displacement vector 

According to COROLLARY 1 proved in (Buchert & Ehlers 1993), a local form can be obtained for second 
order displacements. It reads 

Vo5(2) = Vo5(Ao5) - (Vo5 • Vo) Vo5 ; Vo5 x AoVo5 = . {A.2a,h,c,d) 

The local form (A. 2) is constructed such that its divergence agrees with the source term in (2g), its curl is, 

however, in general non-zero, it only vanishes if (A.2b,c,d) are statisfied. Inserting the potential (A. la) we 

immediately obtain the second-order displacement vector: 

2 (ax sin(27rA')(aj, cos(27rY) + cos(27r2')) \ 
Vo5(2) = — aj,sin(27rY)(a^cos(27rA:) + a^cos(27rZ)) . {A.2e) 
\az SHi{2TTZ){ax cos{2ttX) + ay cos(27rY)) j 



The vector (A.2e) is curl-free as can be easily demonstrated, so it obeys the constraints (A.2b,c,d) necessary 
to admit a potential. This potential can now be guessed from (A.2e) to be of the form 

iS*-^-* := — (ax<^y cos(27rX) cos(27rY) + aj,a^ cos(27rY) cos(27r2') 

+a^a^ cos(27rX)cos(27rZ)) . (-4.2/) 

A. III. The third order displacement vector of the "truncated model" 

Similarily, we can ask for a local vector form whose divergence agrees with the source term in equation (2h). 
An expression given in Buchert (1994, COROLLARY 1) has the required property: The vector Vo5(^") with 
the components 

(Vo5(3")). = ^(Vo5(i)),,J,?, iA.3a) 



has the property 

Ao^^^") = 3777(5J2) , 

where jf^. are the subdeterminants of the tensor (S'''pj,) (a comma always denotes partial derivative with 
respect to Lagrangian coordinates). The following constraints have to be satisfied in order that Voi5'-^"-' be 
curl-free: 

^(Vo5W),4,_^.] = , k^j . {AM,c,d) 

i 

Inserting the potential (A. la) into (A. 3a) gives for the displacement vector 

^3 /sin(27rX)cos(27rY)cos(27rZ)\ 
Vo^^^") = |- a^aj,a^ sin(27rY) cos(27rX)cos(27rZ) . (A.3e) 
\ sin(27rZ) cos(27rX) cos(27rY) / 

Again, the vector (A.3e) is found to be curl-free which renders the contraints (A.3b,c,d) satisfied. A potential 
generating this displacement is again easily found from (A.3e). It reads 

:= -- — - a^aya, (cos(27rX) cos(27rY) cos(27rZ)) . (^-3/) 
(27r)^ 

A. IV. The third order displacement vector of the interaction term 
longitudinal part 

The source term in (2i) which describes the longitudinal part of the interaction of first- and second-order 
perturbations has a similar structure as the second-order source term (2g). We are able to construct a local 
form by analogy (Buchert k Ehlers 1993, COROLLARY 1): 

Ai (Vo5(2)(Ao5(i)) - (Vo5(2) • Vo)Vo5(i)) + A2 {v qS^^X^qS^"^^ - VqS^^^ ■ Vo)Vo5(2)) . (AAa) 

Here, the linear combination of the two possible integrals as a general integral has to be taken, where 
Ai + A2 = 1. In order to satisfy the requirement that the vector (A. 4a) be a solution of the Poisson equation 
(2i), we have to assure that it is curl-free which implies (BE93, COROLLARY 1): 

Ai (Vo5(2) X AoVo5(^)) + A2 (Vo^^^^ x AqVq^^^)) = . {AAb,c,d) 



As can be seen from (A4), we have to determine the parameters Ai and A2 suitably in order to fulfil the 
constraints (A.4b,c,d). Although, we can find the two first integrals for the potential (A. la), the resulting 
vectors are not curl-free. It is a matter of some algebra until one finds the correct linear combination of the 
two vectors, which is curl-free. This can be achieved by first guessing the form of the potential 5'-^'-' from 
the two vectors. It is clear that, in general, we will not be successful. We obtain Ai = |, A2 = |. Thus, the 
displacement vector reads 

3 % ( '^'^ sin(27rA'){aj, cos(27rY) + cos(27r2')]^ \ ^ 
Vo5(^') = |^{- a,sin(27rY)[a,cos(27rX) + a,cos(27rZ)]2 + 
y sm(2'KZ)[ax cos(27rA') + ay cos(27rY)]^ J 

(ax sin(27rA')[aj; cos(27rA')(aj, cos(27rY) + cos(27rZ)) — [ay cos(27rY) — cos(27rZ)]^ + [a'l + a^)] \ 
ay sin(27rY)[aj, cos(27rY)(aj; cos(27rA') + a^ cosi^irZ)) — [a^ cos(27rX) — a^ cos(27r2')]^ + {aj. + a'j,)] j j> 
a^ sin(27r2')[a^ cos(27r2')(aj; cos(27rX) + ay cos(27rY)) — [a^ cos(27rX) — ay cos(27rY)]^ + (a^ + a^)] / 

a^ sin(27rX)[2aj; cos(27rX)(aj, cos(27rY) + a^ cos(27r2')) + IQaya^ cos(27rY) cos(27r2')] \ 
ay sin(27rY)[2aj, cos(27rY)(aj; cos(27rX) + a^ cos(27r2')) + IQa^a^ cos(27rX) cos(27r2')] 
a^ sin(27r2')[2a^ cos(27r2')(aj; cos(27rX) + ay cos(27rY)) + IQa^ay cos(27rX) cos(27rY)] / 

/ a^ sin(27rX)[a2 cos2(27rY) + al cos2(27rZ) + 2{a'l + a^)] \ 
+ ttj, sin(27rY)[af cos2(27rX) + cos2(27rZ) + 2(af + a^)] | , {AAe) 
\ a, sin(27rZ)[a2 cos2(27rX) + a'^ cos2(27rY) + 2{al + a^)] / 

with the potential 

:= --j^-y g [a't cos2(27rX)[aj, cos(27rY) + a, cos(27rZ)] + IQa^aya, cos(27rX) cos(27rY) cos(27rZ) 

-\-a'y cos^(27rY)[aj; cos(27rX) + cos(27r2')] + a'j, cos^(27r2')[aj; cos(27rX) + ay cos(27rY)] 

+2[ax{al + al) cos(27rX) + ay{al + a^) cos(27rY) + a^{al + a^) cos(27rZ)]| . (AAf) 




27r5 



A.V. The third order displacement vector of the interaction term 
transverse part 

Finally, we ask for a first integral of the transverse part of the interaction vector (2j,k,l). In (Buchert 1994, 
COROLLARY 2) the vector form needed has been given again as a linear combination of the two possible 
integrals. The vector 

S := -Vo X 5(^<=) = 

111 ((Vo5(2) . Vo)Vo5(i)) - 112 ((Vo5(i) • Vo)Vo5(2)) (A.5a, h, c) 

has the property 

d{S^^.\S^;p\x,) 
(Vo X S)fc = fp^y-^^^-^-^-^ , i,i,k= 1,2,?, {cyclic) . 

We have to assure jii + = 1. In order to satisfy the requirement that the vector components (A.5a,b,c) be 
solutions of the Poisson equations (2j,k,l), we have to guarantee that the vector field S is source-free which 
implies after using well-known vector identities 



Ao(Vo5(i)Vo5(2))+ 



Vo5(2)AoVo5(i) - Vo5(i)AoVo5(2)) +//2 (Vo5(2)AoVo5(i) - Vo5(i)AoVo5(2)) = . (A.bd) 



Again, we find tlie two integrals after inserting the potential (A. la) to be not source-free. We have to 
determine the correct linear combination of the two integrals. As in the longitudinal case we first guess the 
form of the vector potential 5'-^'^-' from the two integrals. We are successful with the parameters /^i = |, 
and = f, but we have to add another function jipV oV^^^; jjip = such that the total displacement is 
source-free. The potential V^^^ is given by 

•p(i) = a^ia'y + a^) cos(27rX) + ay{al + a]) cos(27rY) + a,{al + a'y) cos(27rZ) . (A.5e) 

(This is possible, since the local forms discussed above are only determined up to the gradient of some 
potential, see Buchert 1994). The vector displacement S = — Vo x S^^""^ reads 

3 g / sin(27rX) cos(27rX)[aj, cos(27rY) + cos(27r2')] \ ^ 
S = |-|- a2sin(27rY)cos(27rY)[a^cos(27rX) + cos(27rZ)] - -• 
^ ^ \al sin(27rZ) cos(27rZ)[a^ cos(27rX) + ay cos(27rY)] / ^ 

(al sin(27rX)[cos(27rX)(aj, cos(27rY) + cos(27rZ)) + a| cos2(27rY) + aj, cos^CIttZ) - {a 
a'l sin(27rY)[cos(27rY)(a^ cos(27rX) + a, cos(27rZ)) + af cos2(27rX) + al cos'^{2ttZ) - {a 
al sin(27rZ)[cos(27rZ)(aj; cos(27rX) + ay cos(27rY)) + a]. cos2(27rX) + cos2(27rY) - (a 

^3 Y +a?)sin(27rX)\ ^3 ^ 

- TTI 7 %(4 + a?)sin(27rY) = ^ T' 
\ a,^(a^ + afj) sm(27rZ) / 

/ a'l sin(27rX)[cos(27rX)(aj, cos(27rY) + a, cos(27rZ)) - 2a'^ cos2(27rY) - 2al cos2(27rZ) + (a^ + aj)] ' 
a'l sin(27rY)[cos(27rY)(a^ cos(27rX) + a, cos(2ttZ)) - 2al cos^(2ttX) - 2al cos'^{2ttZ) + (a| + al)] 
\ al sin(27rZ)[cos(27rZ)(a^ cos(27rX) + ay cos(27rY)) - 2al cos'^{2ttX) - 2al cos'^{2ttY) + {al + a^)] ^ 




(^•5/) 



with the vector-potential 



3 I / ttj/tt^ sin(27rY) sin(27r2')(aj, cos(27rY) — cos(27r2')) \ 
5(^<=) := - a^a^ sin(27rX)sin(27rZ)(a^cos(27rZ) -a^cos(27rX)) . (A.bg) 
^ ' \axay sm(2TrX) sm(2TrY)(ax cos(2TrX) — ay cos(2TrY)) J 



A. VI. Remarks 

For the following discussion we need the explicit expressions of the source terms in the solution (2) for the 
special model (A. la). We derive 

7(5 i_fc) = ejaj; cos(27rX) + ttj, cos(27rY) + cos(27r2')} , (A. 6a) 

Il(S^i^k) = s^io^xCty cos(27rX) cos(27rY) + aya^ cos(27rY) cos(27r2') + a^a^ cos(27rX) cos(27r2')} , (A.6b) 
III{S^i^k) = e^{a^aya, cos(27rX) cos(27rY) cos(27rZ)} , (A. 6c) 

E = ^"i"' cos2(27rX)K cos(27rY) + a, cos{2^Z)] 

a,b,c 

-\-a'y cos^(27rY)[aj; cos(27rX) + cos(27r2')] + al cos^(27r2')[aj; cos(27rX) + ay cos(27rY)] 

+6aj;aj,a^ cos(27rX) cos(27rY) cos(27r2')} , {A.Qd) 



d{S^?,S,p,Xq) ( sm(27rY) sm(27rZ)(aj, cos(27rY) - a, cos(27rZ)) \ 

'^pq\j-^7^ — Y Y \ ~ \ "2^"^sin(27rX)sin(27rZ)(a^ cos(27rZ) - aj;Cos(27rX)) . {A.&e,f,g) 
C(Ai , A2, Aaj \^ ^^^^ sin(27rX) sin(27rY)(a^ cos(27rX) - ay cos(27rY)) / 

From the generating functions constructed above we infer the following property: except for the longitudinal 
part of the interaction term, the perturbation potentials obey equations which are typical for bound systems: 

Ao5(i) = -(27r)2 5(1) , {A.la) 

Ao5(2) = -(27r)2 2 S'-'^'^ , (A.lh) 

AoS^^"^ = -(27r)2 3 S^^"^ , (A. 7c) 

Ao5(3') = -(27r)2 {5 5(3') - 4 5(3") + 2 •p(i)} , (A.ld) 

Ao5(3<=) = -(27r)2 5 S'^^'^ . (A.le) 



Recall that the condition (A. 7a) implies Voi5(to) cx u for an initially irrotational peculiar-velocity field u(to)', 
Vo X u(to) = 0, i.e., the motion is initiated to follow the gradient of the density-contrast field. At the third 
order the evolution model shows that this property of the fiow is lost. 

(The algebraic program we used to compute the perturbation potentials for Model II also reproduces the 
perturbation potentials derived here.) 



